The shape of a moving fluxon in stacked Josephson junctions. 
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We study numerically and analytically the shape of a single fluxon moving in a double stacked 
Josephson junctions (SJJ's) for various junction parameters. We show that the fluxon in a double 
SJJ's consists of two components, which are characterized by different Swihart velocities and Joseph- 
son penetration depths. The weight coefficients of the two components depend on the parameters 
of the junctions and the velocity of the fluxon. It is shown that the fluxon in SJJ's may have an un- 
usual shape with an inverted magnetic field in the second junction when the velocity of the fluxon is 
approaching the lower Swihart velocity. Finally, we study the influence of fluxon shape on flux-flow 
current-voltage characteristics and analyze the spectrum of Cherenkov radiation for fluxon velocity 
above the lower Swihart velocity. Analytic expression for the wavelength of Cherenkov radiation is 
derived. 
I. INTRODUCTION 



Properties of stacked Josephson junctions (SJJ's) are of considerable interest both for applications in cryoelectronics 
^2 , an d for fundamental physics. A particular interest in SJJ's was stimulated by the discovery of high-T c superconductors 
(HTSC). Highly anisotropic HTSC compounds, such as k^S^CaC^Os+i, may be considered as stacks of atomic 



scale intrinsic Josephson junctions [[l). The layered structure determines many of the unusual properties of HTSC. 
Behavior of model low-T c SJJ's and HTSC exhibit many similarities Due to mutual coupling of the junctions 
in the stack, the physical properties of SJJ's can be qualitatively different from those of single Josephson junctions 
(JJ's). Therefore, a one-to-one comparison between single and stacked Josephson junctions is difficult to do. Hence, 
the basic properties of SJJ's have to be studied in order to describe correctly the Josephson behavior of layered 
O ' superconductors. 

Perpendicular (c-axis) transport measurements in magnetic field, H, parallel to layers (afr-plane) is an explicit 
way of studying Josephson phenomena in SJJ's. In this case the magnetic field penetrates the stack in the form of 
t-H , Josephson-type vortices (fluxons), and the c-axis voltage is caused by motion of such vortices along the layers. The 
shape of a fluxon in SJJ's is different both from that of an Abricosov vortex in bulk superconductor, since it does not 
have a normal core, and from that of single J J, since the circulating currents are not confined within one junction. The 
behavior of SJJ's becomes particularly complicated when the length of the stack in one direction, L, is much longer 
, than the Josephson penetration depth, Aj. One of the unusual properties of long SJJ's is the existence of multiple 
quasi-equilibrium fluxon modes || , and submodes j|] , which are characterized by different fluxon configurations in the 
stack. Due to the existence of such modes/submodcs, the state of the stack is not well defined by external conditions. 
Rather it can be described only statistically with a certain probability of being in any of the quasi-equilibrium states. 
Experimental evidences for the existence of such modes were obtained both for HTSC intrinsic SJJ's j|-[6) and low-T c 
multilayers [Q. The existence of fluxon modes/submodes dramatically changes the behavior of long strongly coupled 
SJJ's with respect to that of single long JJ's. An example of this is the critical current, I c , which becomes multiple 
valued [Q-)6| ; the fluctuations of I c become anomalously large ||,[| ; and the magnetic field dependence of I c becomes 
C , very complicated without periodicity in H |Q| . 

For understanding both the static and dynamic properties of SJJ's, the shape of the fluxon in SJJ's is important, 
and should be determined. In the static case, the shape of the single fluxon was studied for layered superconductors 
consisting of an infinite number of thin identical Q or nonidentical || layers and for SJJ's PUlPf]. In our previous work 
||, we have shown that in double SJJ's, two special single component fluxon solutions exist, which are characterized by 
different Swihart velocities and Josephson penetration depths. An approximate analytic fluxon solution was suggested 
as a linear combination of the single component solutions. For the static case, the approximate solution was shown 
to be in a quantitative agreement with numerically obtained solutions. Extending the approximate analytic solution 
to the dynamic case, it was shown that drastic changes in the fluxon shape could occur with increasing the fluxon 
velocity, resulting e.g. in possible inversion of the sign of the magnetic field in the second junction and appearance of 
attractive fluxon interaction |3|. On the other hand, the choice between the special single component solutions and 
the approximate analytic fluxon solution was not addressed and the dependence of the fluxon shape on the junction 
parameters was not studied. Using the perturbation approach, the second order correction to the approximate analytic 
solution was derived and the accuracy of the solution was analyzed in the recent paper |lJj] . 

To our knowledge, no comprehensive analysis of the single fluxon shape in SJJ's exists for the dynamic case. The 
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scope of the current paper is to study quantitatively the shape of the moving fluxon in double SJJ's for various 
junction parameters. Our analysis is based on numerical simulations and analytical treatments of the coupled sine- 
Gordon equation, which describes the physical properties of SJJ's PJ. We show that the single moving fluxon in 
double SJJ's may be described by both a single component solution and a double component solution depending on 
the parameters of the stack and the fluxon velocity. Moreover, the shape of the fluxon may be quite anomalous, with 
inverted magnetic field in the second junction and with nonmonotonous change of phase. 

The paper is organized as follows: in section II, we rewrite the coupled sine-Gordon equation for the case of 
solitonic fluxon motion and review analytic single fluxon solutions obtained in Refs. PJH]|. In section III, we present 
the results of numerical simulations for frictionless fluxon motion for different parameters of SJJ's, compare it with 
analytical predictions of Ref. we also formulate and verify conditions for observation of different fluxon solutions. 
In section IV, we discuss implementations of the fluxon shape in experimental situation. In subsection IV. A, we study 
the influence of finite damping and simulate current-voltage characteristics. Finally, in subsection IV. B we consider 
the case of nonsolitonic fluxon motion with the propagation velocity larger than the lower Swihart velocity. We have 
shown that such fluxon motion is accompanied by plasma wave exitations and derive the expression for the wavelength 
of such " Cherenkov" radiation. 

II. GENERAL RELATIONS 

We consider a double stack with the overlap geometry, consisting of junctions 1 and 2 with the following parameters: 
J C i -the critical current density, Ci -the capacitance, ti- the thickness of the tunnel barrier between the layers, di and 
Xsi - the thickness and London penetration depth of superconducting layers and L- the length of the stack. Hereafter 
the subscript i on a quantity represents its number. The elements of the stack are numerated from the bottom to the 
top, so that junction i consists of superconducting layers i, i + 1 and the tunnel barrier i. The fluxon will be placed 
in junction 1, if not stated otherwise. 

The physical properties of SJJ's are described by the coupled sine-Gordon equation pJ, which for a double stack 
with the overlap geometry can be written as: 
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where tpi j2 are gauge invariant phase differences in JJ's 1 and 2, 'prime' and 'dot' on the quantity represent 
partial derivatives in space and time, respectively. Space and time are normalized to Josephson penetration depth, 

J 1 = V 8ttV c iAi ' and inverted Josephson plasma frequency, u^ 1 = J ■ 
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respectively, of the single J J 1. Here 



$0 is the flux quantum, c is the velocity of light in vacuum and 
A, = U + Xsicoth (^-) + Xsi+icoth (^7) , 

Si = X Si cosech (j^J ■ 

The last terms in the right hand side of Eq.(l) represent total currents in the JJ's, which consist of superconduct- 
ing, displacement and quasiparticle contributions, and Jb represents bias current density. Viscous damping due to 

— 1/2 

quasiparticle current is characterized by the damping coefficient, a, = p ci , where p c i is the McCumber parameter 
of single J J i. 

The coupling strength in the double SJJ's is described by a coupling parameter 

S = ^ 2 
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where Hq 

For the soliton-like fluxon motion, the phase differences in the stack remain unchanged in the coordinate frame 
moving along with the fluxon. Introducing the self-coordinate of the fluxon, £ = x — ut, and neglecting damping 
coefficient, we simplify Eq.(l) and rewrite it as a system of coupled ordinary differential equations (ODE): 
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and coi = AjjWpi is the Swihart velocity of the single junction 1. Comparing Eqs. (1) and (3) it is seen that solution 
of the coupled sine-Gordon equation for soliton-like fluxon motion is now reduced to solution of the static problem, 
but with parameters depending on the fluxon velocity. 
Eq.(3) has a first integral, 
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which reduces to that from Ref. j3j for the static case, u = 0. Here C is a constant of the first integral. 
A. Special single component solutions 

In Ref . fel it was shown that Eq.(l), linearized with respect to ip 2 , allows two special single component solutions of 
the type fjj] 



<Pi(0 = F i,2 = iarctan [exp (£/Ai, 2 )] 
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where ki j2 are solutions of the quadratic equation: 
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Therefore, for a double SJJ's there exist two characteristic Josephson penetration depths, 
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and two characteristic velocities 
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_B. Double component solution 

Taking the single component solutions as eigen-functions of the linearized coupled sine-Gordon equation, an ap- 
proximate analytic single fluxon solution in J J 1 was obtained in Ref. || : 
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Here are the single component solutions, Eq.(6). Recently this solution was rederived more rigorously in Ref. 
jnj. It was shown, that for the static case Eq.(10) gives perfect approximation for ipi in the whole space region and 
for arbitrary parameters of the stack |11 M. Using the perturbation approach the second order correction to Eq.(10) 
was obtained in Ref. |11|. As it is seen from Eq.(10) the single fluxon in double SJJ's consists of two components. 
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From Eq. (8) it is seen that both components contract with increasing velocity, but the characteristic velocities for the 
contraction are different for each component and are given by Eq.(9). For identical junctions the contraction of each 
component is of the Lorentz type, however, the contraction of the fluxon itself is different from Lorentz contraction. 
This is a consequence of the absence of Lorentz invariance for the coupled sine-Gordon equation. For nonidentical 
junctions, the parameters ki,2, depend on the fluxon velocity and thus contraction of the components is somewhat 
different from Lorentz contraction. In this case the maximum characteristic velocity should be obtained from the 
equation u — c\ 2 and Eqs. (7,9). By analogy with single JJ's we'll refer the maximum characteristic velocities to as 
Swihart velocities, c\^- In the general case, Swihart velocities are equal to p3[ 

\/2coiCo2 

Cl,2 = < (11) 



\j 4i + c§ 2 ± \f (cqi ~ c 2 Q2 f + 45^P 



where C02 = c oiy cIaI * s ^ e Swihart velocity of the single JJ2. 

The most crucial changes in the fluxon shape occur as the velocity approaches the lowest Swihart velocity, u — > c± . 
Then the first component is totally squeezed, Ai — * 0, while contraction of the second component remains marginal, 
see Eq.(8). In this case the two components become clearly distinguishable: the F\ component transforms into a 
step-like function which changes from zero to 2tt within the distance Ai at the fluxon center, while outside the central 
region the shape of the fluxon is defined by the F 2 component. From Eq.(10) it follows that 

sin^px) [k 2 ,\x\ > Ai; 

(u -> ci) . (12) 



sin(tp2) \—Kx,x—*0 
For q 2 ^ 2 — 1, and u = e"i, the parameters Ki t 2 are equal to: 




The parameters K\ : 2 determine the weight coefficients of the components. From Eqs. (10, 13) it follows that F\ com- 
ponent dominates in junction 1 for J c2 j ' J c \ <C 1, and Fq, dominates for J&j J c i 3> 1. 

From the analysis above it is seen that the linearized coupled sine-Gordon equation allows both the single component 
solutions, Fi.2, Eq.(6) and the double component solution, Eq.(10). At this stage it is not clear which of the solutions, 
Eqs. (6, 10), should be realized in SJJ's, since all three solutions have the same accuracy with respect to Eq.(l). In 
Refs. H[[l] it was shown that it is the double component solution Eq.(10) which is realized in the static case. However, 
it was suggested that a single component solution could be achieved at high fluxon velocities. Indeed, as we will show 
below, in the dynamic case both single and double component solutions can exist and even coexist, depending on 
parameters of the stack and the fluxon velocity. What is important, however, is that these are always the components 
Fi.2 described by Eqs. (6,7) which constitute the fluxon. 

III. FRICTIONLESS CASE 

In this section we will consider unperturbed, on = Jt, = 0, frictionless fluxon motion. We analyze the pure solitonic 
fluxon motion for various junction parameters, make general conclusions about transformation of the fluxon shape 
in dynamics and compare it with analytical predictions. The effect of finite damping and bias will be considered 
in section IV. For numerical simulations of frictionless fluxon motion we used Eq.(3). The numerical procedure was 
based on a finite difference method with successive iterations of ODE in Eq.(3). The boundary conditions were such 
that the total phase shift is equal to 2?t in the junction containing a fluxon and zero in the other one. The fluxon 
will be placed in JJ1 if not stated otherwise. We will consider five cases: A) a stack of identical junctions, B) stack 
with different critical current densities, C) fluxon in the junction with lower critical current density, D) fluxon in the 
junction with the higher critical current density, E) a stack with difference both in critical currents and electrodes. 
Finally in subsection III.F we derive and verify conditions for observation of single and double component fluxon at 
u — > c i for various parameters of SJJ's. 

A. Identical junctions: double component solution. 

In Fig.l profiles of a) phase differences <£>i,2j b) the ratio sin(ipi) / sin(ip2) , and c) magnetic inductions B\ 2 of a 
single fluxon in junction 1 are shown for a double stack consisting of identical strongly coupled junctions and for 
different fluxon velocities, u/cx=0, 0.61, 0.92, 0.98, 0.998, 0.9999 (from left to right curve). The curves were shifted 
for clarity along the :r-axis. Parameters of the stack are: di =U = O.OIAji, Xsi = O.IAji, S ~ 0.5, C% = Ci- In Fig.l 
a) dotted lines show profiles obtained from the analytic double component solution Eq.(10). Solid and dashed lines 
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in Fig. 1 a,c) represent results of numerical simulation for junctions 1 and 2, respectively. The data in Figs.l b,c) are 
obtained numerically. The magnetic induction is normalized to Hq = ■ 

As the velocity approaches the lower Swihart velocity, ci, the existence of the two fluxon components becomes 
clearly seen. For identical junctions, as it follows from Eqs.(10,13), exactly one half of the fluxon belongs to each 
component. The F\ component transforms to a one-7r step. Outside the fluxon center the fluxon is defined entirely by 
the F2 component, which is only marginally contracted. Moreover, from Fig.l a) it is seen that at u ~ c\ the phase 
differences in both junctions are equal outside the fluxon center in agreement with analytical prediction, Eq.(12). This 
is illustrated in Fig. 1 b), from which it is seen that the ratio sm(^) / sin(if2) approaches unity as it — > ci. From 
Fig.l a) it is seen that the approximate analytic solution is in good agreement with numerical solution for all fluxon 
velocities. 

Another unusual feature of the moving fluxon in SJJ's is seen from Fig. 1 c). A dip in E>2 is developed with 
increasing fluxon velocity, leading to inversion of the sign at high velocities. Such behavior was predicted analytically 
in Ref . |3|] . Here we confirm the existence of this phenomenon by numerical simulation. 

B. Nonidentical junctions: transformation from F 2 to i*\ component solution with decreasing J c il^c\ : 

In Ref. H it was shown that it is the double component solution that is realized in the static case for arbitrary 
parameters of the stack. From the numerical analysis of the dynamic case we have found that the fluxon shape is well 
described by the double component solution up to velocities very close to 7i\. This is illustrated by Fig. 2, in which the 
shape of the fluxon in junction 1 moving with the velocity close to the lower Swihart velocity, u — 0.98ci, is shown 
for different critical currents, J C 2 / Jci — ^0', 2; 1; 0.5; 0.1, from left to right curve. The rest of the parameters of the 
stack and the way of presentation is the same as in Fig. la). From Fig. 2 it is seen that the fluxon shape is in a good 
agreement with the analytical double component solution, Eq.(10), up to velocities very close to c%. Note, that the 
agreement is even better at lower fluxon velocities. From Fig. 2 it can be seen that a gradual transformation of the 
fluxon shape from the uncontracted F2 component solution to the contracted Fi component solution takes place as 
Jci/Jci decreases. In terms of the analytical double component solution, Eq.(10), this is caused by a decrease of the 
weight coefficient K2 , Eq.(13), of the F2 component. 

As the velocity approaches the lower Swihart velocity, ci, the fluxon shape can remain the double-component, as 
shown in Fig.l, or a transformation of the fluxon shape could take place. As we have found from our numerical 
simulations the transformation and the final shape of the fluxon at u = c\ strongly depend on parameters of the 
stack. 

C. Fluxon in a weaker junction: Uncontracted single component solution. 

In Fig. 3 the fluxon shape is shown for the case when the fluxon is placed in the weaker junction, J C 2/ Jci — 2, the 
rest of the parameters and the way of presentation are the same as in Fig.l. As it is seen from Fig. 3 at velocities up 
to 0.98ci the shape of the fluxon is well described by the double component solution, Eq.(10). At higher velocities 
transformation to the single F2 component solution, Eq. 6, with takes place. This is illustrated in Fig. 3 b), from 
which it is seen that as u — » ci, sin(ipi) j sin(if2) - *• 1^2 — 2 in the whole space region. From Fig. 3 c) it is seen that a 
dip in E>2 is reduced with respect to that in Fig.l c), due to absence of the squeezed F\ component. 

D. Fluxon in a stronger junction: two component solution. 

In Fig. 4 the fluxon shape is shown for the case when the fluxon is placed in the stronger junction, J c2 / J c \ = 0.5, the 
rest of the parameters and the way of presentation are the same as in Fig. 1. It is seen that at velocities up to 0.98ci 
the shape of the fluxon is well described by the double component solution, Eq.(10). At higher velocities the fluxon still 
has contracted and uncontracted components -F1.2. The existence of the two fluxon components is clearly seen from 
Fig. 4 b). As u — > ci, sin(ipi) / sin(ip2) K2 — 0.5 outside the center of the fluxon and sin(ipi) / sin((p%) — * — K\ = 1 
in the center, in agreement with Eqs.(12,13). However, transformation of the fluxon shape with respect to Eq.(10) 
takes place. From Fig. 4 a) it is seen that in the left half-space (i.e. from far left to the immediate left of the fluxon 
center) the phase shift in junctions 1 and 2 approaches zero and it, respectively, and belongs to the uncontracted F2 
component, as seen from Fig. 4b). Therefore, there is a single F 2 component fluxon placed in junction 2 (the weaker 
junction). The situation in the left half space is then analogous to that in Fig. 3. This is shown in Fig. 5, in which we 
replotted the curves at u — 0.9999ci from Figs. 3a) and 4a). Solid and dashed curves in Fig. 5 represent ipi and p 2 , 
respectively, for the fluxon in the stronger junction from Fig. 4, while dashed-dotted and dotted curves represent tpi 
and (f2, respectively, for the fluxon in the weaker junction from Fig. 3. Note, that in Fig. 5 we normalized the x-axis 
to Aji for the case J C 2/J c i = 2 and to Xj2 for the case J C 2/Jci = 0.5, since the F% component is then situated in 
junction 2. From Fig. 5 it is seen that after such rescaling (^2,1 from Fig. 4a) merge with tpi^ from Fig. 3a) in the 
left half-space. In the central region a step-like change of phase shift takes place in both junctions. In junction 1 
the phase jumps on +2-7T, which means that there is a single Fi component fluxon and in junction 2 the phase drops 
on — 2-7T, representing the single F\ component antifluxon. The overall fluxon shape at u = ci for the case when the 
fluxon is in the stronger junction, Fig. 4, can be written as: 
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<p 1 =F 1 + Image((p 2 ) ^ 
¥>2 = F 2 - F 1 

so that in junction 1 there is the contracted single component F\ fluxon plus an image in a form of a ripple from 
junction 2 and in junction 2 there is an uncontracted fluxon F 2 - contracted antifluxon F\ pair. The total phase shifts 
in junction 1 and 2 are 27r and zero, respectively, however, the phase shift in junction 1 increases nonmonotonously 
and has two local maxima and minima, see Fig. 4 a). From Fig. 4 c) it is seen that the dip in B 2 in this case is even 
more pronounced than that for identical SJJ's, Fig.l c). This is due to the increase of the weight coefficient of the 
squeezed F\ component. 

Fig. 6 shows profiles of the fluxon moving with velocity very close to the lowest Swihart velocity, u = 0.9999ci, for 
different critical current densities J c2 / 'J c i = 10-7- 0.1 increasing sequentially from the left to the right curve. The rest 
of the stack parameters and the way of presentation is the same as in Fig. 1. From Fig. 6 it is seen how the shape of 
the fluxon is changed with J c2 /J c i- When the fluxon is placed in the weaker junction, the fluxon shape at the lowest 
Swihart velocity is described by the single F 2 component. For the case of Figs. 1-6, C 2 /C\ = A 2 /Ai = 1, so that n 2 = 
■JY_ In Fig. 6 b) the dependence sin(tpi) / sin((p 2 ) — -j^ is clearly visible in the whole space region for J c2 /J c i > 1. 
When the fluxon is placed in the stronger junction, it has two components, F\ t2 . From Fig. 6 b) it is seen that for 
the case J c2 /J c i < 1, the fluxon shape in the center is determined by the F\ component, sin(ipi) / sin(ip 2 ) = —K\ = 1, 
while outside the center the shape is given by the F 2 component, sinUpi) / sinUp 2 ) — k 2 — 4 s2 -, confirming that the 
parameters of the single component solutions are always those that predicted by Eqs.(6,7). 

From Fig. 6 it is seen that the transition from a double to a single component solution for J c2 /J c i > 1 is gradual. 
Outside the fluxon center this transition is well described by a gradual increase of the weigh coefficient of the F 2 
component, Eq.(10,13). In the center the exact shape of the fluxon can be obtained from the first integral, Eq.(5). 
For the case = 1, at u = C\ the first integral reduces to: 



J c iA lC os(v»i)\ 2 (vi*) J c2 
1 + ^— r -, — - ' = 1 - cos + — (1 - cos (<p 2 )) , (15) 



1 - S 2 V ^A 2 cos (ip 2 ) J 2 Jcl 

for the F 2 single component solution. From Eq.(15) it is seen that at the fluxon center, x = 0, the effective Josephson 
depth is equal to: 



-M'-SeJVi^ (16) 

so that X e ff gradually increases from zero to Xji^J jA^z,aJS becomes larger than unity. The inequality 

Jc2A 2 



JciAi 



> 1, (17) 



is then a necessary (but as we'll show below not sufficient) condition for the existence of the single component F 2 
solution at u = c±, for = 1, since then \ e ff remains finite at u — C\. 

On the other hand, the transition from a double component solution, Eq.(10), to the solution, Eq.(14), for J c2 /J c i < 
1 is sharp. However, the closer is J c2 /J c \ to unity, the closer must be the fluxon velocity to C\ in order to observe 
this transformation, as it can be seen from Figs. 4,6. 

E. Nonidentical electrodes: Bifurcations and more complicated two component solutions. 

So far we have considered the case when only critical current densities of JJ's in the stack were different. Another 
common type of the difference between JJ's is the difference in the properties of electrodes. In Fig. 7 the fluxon shape 
is shown for the case, A 2 ~ 2.5Ai. Physically this means that the third electrode has either larger London penetration 
depth, Ag3 = 2Aji j2 , or it is thinner than the rest of the electrodes, d\_ 2 — Ad%, see definitions in sec. II. The rest of 
the parameters are J c2 /J c i = 0.5, d\. 2 = U = O.OIAji, Xsi,2 = O.IAji, S ~ 0.31, = i. The way of presentation 

is the same as in Fig. 1. As it is seen from Fig. 7, at velocities up to 0.98ci the shape of the fluxon is well described 
by the analytic double component solution, Eq.(10). From Fig. 7b) it is seen that outside the fluxon center the phase 
distribution is determined by the F 2 component with k 2 ~ 0.79 given by Eq.(13). However, at u <~ 0.998ci the system 
exhibit bifurcations and a sudden switching to the solution given by Eq.(14) occurs. At slightly larger velocity another 
bifurcation takes place resulting in switching to yet another solution. 

The switching between the solutions is hysteretic. If we start reducing the fluxon velocity, the switching back to 
the double component solution takes place at somewhat lower velocity. Therefore there is a certain region of fluxon 
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velocities for which those solutions coexist. Such situation is illustrated in Figs. 8 a-c), in which phase distributions 
<pi.2 for three possible solutions are shown for u ~ 0.998ci and for the same stack as in Fig. 7. As usual, solid and 
dashed curves in Fig. 8 represent ipi and f 2} respectively, obtained from numerical simulations. 

Fig 8 a) shows the solution coming from low velocities. Dashed-dotted and dotted curves in Fig. 8 a) represent ip\ 
and ip2, respectively, given by the approximate double component solution, Eq.(lO), showing good agreement with 
numerical solution. 

With increasing velocity, switching to the solution of the type given by Eq.(14) takes place. This solution is shown 
in Fig. 8 b). This solution exist only in a very narrow velocity interval and with further increase of velocity a switching 
to a new solution occurs. The new solution is shown in Fig. 8 c). It is characterized by a large phase variation in the 
second junction. Such solution exists up to ci. 

Looking at the fluxon shape at u = 0.9999ci from Fig. 7 we see that it consists of three parts: (i) at the fluxon 
center, x — 0, in junction 1 there is a 2-k step, which as we will show in a moment, belongs to the pure F\ component 
and in junction 2 there is a At: (!) drop belonging to two flux quantum, — 2iq, antifluxon. (ii) At x — xq ~ 0.7Aji 
and (hi) x = —xq in junction 2 there is a 2ir increase and in junction 1 there is a ripple like phase shift. As it can be 
seen from Figs. 7 a,c) the features at x = ±x contain both contracted and uncontracted parts and are in fact given 
by the double component solutions, Eq.(10). This is illustrated in Fig. 9 a) in which solid and dashed lines show phase 
distributions ip± and ip 2 , respectively, from Fig. 7 a) for u = 0.9999ci and dashed-dotted and dotted curves represent 
ipi and <f2, respectively, given by the analytic double component solution, Eq.(10), shifted by — x along the x- axis. 
Obviously, the features at x = ±xq correspond to the double component soliton placed in junction 2. Due to spatial 
separation between the contracted centrum of the fluxon and the double component features at x = ±xq, we can 
analyze the shape of the central contracted part. From Fig. 7 c) it is seen that the magnetic induction in junction 1 
at x = 0, -Bi(O), increases sharply as the velocity approaches c\. If the phase distribution in the central region is given 

by the F\ component, then from Eq. (2) it follows that £?i(0) <~ F{(0) ~ A7 1 ~ ( 1 — % j , representing Lorentz 
contraction of F\ component at the lowest Swihart velocity, c\. Fig. 9 b) shows the inverse value of -Bi(O) versus 

( A 1/2 

the Lorentz factor 1 — % ] , for the high velocity solution from Fig. (7). Dots represent the numerically obtained 



values, solid line is the apparent linear fit. Lorentz contraction of the central region at the lowest Swihart velocity, c\, 
is clearly seen from Fig. 9 b), therefore confirming that the central region is given by the pure F\ component. The 
overall shape of the soliton in this case can be described as: 



{ 



if! = Fx + Image{ip 2 ) . . 

<?2 = F 1+2 (x - x ) - 2Fi + F 1+2 (x + x ) ' 1 1 



where F i+2 is the double component solution, given by Eq.(10). From Eq.(18) it is seen that the overall phase shift 
is 2n in junction 1 and zero in junction 2, however the phase growth in junction 1 is nonmonotonous. 

F. Conditions for the existence of single and two component solutions. 

For practical applications of SJJ's in flux-flow oscillators, the shape of the fluxon at the highest propagation velocity 
is crucial. It is then important to know how the shape of the fluxon at u = c\ depends on the parameters of the 
stack. Eq.(17) formulates the necessary condition for the existence of uncontracted single component F 2 solution at 
u = ci, for Q 2 ^ 2 = 1, since then X e ff remains finite at u = C\. However this condition is not sufficient. Indeed, for the 
case of Figs. 7, 8, j^T — 1-25, i- c - Eq.(17) is satisfied. However, we obviously do not observe the single component 
F 2 solution at u = ci, but rather a switching to more complicated two component solutions, Eqs.(14,18), takes place. 
From Fig. 8 a) it is seen that the switching takes place when maximum value of ip 2 approaches 7r/2. If the solution 
were to stay at the special single component F 2 solution, then according to Eq.(6), the maximum value of sin(<p 2 ) 
would be equal to 1/k 2 . This gives us an additional condition for the observation of the single component F 2 solution: 

k 2 > 1, (19) 

In Fig. 10, regions of the existence of the two component (shaded area) and the single component F 2 solutions at 
u = ci are shown for _ \ Numbers in Fig. 10 show the number of components obtained numerically. Solid and 
dashed lines represent the conditions Eq.(19) and Eq.(17), respectively. Arrows indicate the cases considered in Figs. 
1,3,4,7. It is seen that the single component solution exists when both conditions, Eqs.(17,19), are satisfied. 

IV. IMPLEMENTATION FOR EXPERIMENTAL SITUATION 

In the previous section we have studied the unperturbed fluxon motion, on = 0, Jb = 0. This is an idealized case. 
In real experimental situation a* 7^ and 7^ 0. To quantitatively study the influence of damping and bias on 
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the fluxon shape and IVC's, we performed numerical simulation of the coupled sine-Gordon equation, Eq.(l), with 
the dissipation and bias terms. Here we used two different approaches: (i) Considering solitonic-type fluxon motion, 
<fii = we derive the system of one-dimensional ordinary differential equations (ODE), similar to Eq.(3), but with 

dissipation terms, a, 7^ 0, and bias terms, and used the same numerical procedure to solve it. (ii) Alternatively, we 
directly integrated the system of partial differential equations (PDE), Eq. (1) using explicit finite difference method. 
Both approaches have certain advantages and disadvantages. Using ODE, it is possible to calculate fluxon shape for 
arbitrary small damping, while PDE require relatively large damping coefficient. In addition, ODE in comparison 
to PDE, does not have problem with accumulation of error and does not require a long relaxation and averaging 
times. For solitonic fluxon motion both approaches give identical results. On the other hand, ODE are restricted 
to the study of solitonic motion, while the PDE allow more complicated solutions. For simulation of PDE, periodic 
boundary conditions were used, which correspond to fluxon motion in annular SJJ's with L = 10 -j- IOOAji. 
A. Effect of damping and current-voltage characteristics. 

First of all fluxon shape will affect the shape of the flux-flow current-voltage characteristics (IVC's) caused by 
fluxon motion in the stack fill . The shape of flux- flow IVC is determined by a balance between the input power 
to the system from the current bias source and the power dissipated due to a finite damping, see e.g. Ref. In 
a single JJ's as the fluxon velocity approaches the Swihart velocity, Lorentz contraction takes place and the fluxon 
energy increases sharply. The fluxon velocity asymptotically approaches the Swihart velocity with increasing current. 
In the IVC this result in the appearance of an almost vertical step at the velocity matching condition |l(J. Since the 
existence of this step is closely related to the Lorentz contraction of the fluxon, we expect that the step at the velocity 
matching condition should also exist in SJJ's with the characteristic velocity equal to the lower Swihart velocity, c\, 
whenever the fluxon contains the Lorentz contracted part F± . On the other hand, when a pure F2 component solution 
takes place, the fluxon will reach u = ci at a finite current and flux-flow IVC should have a finite slope at u — c±. 
For the case when the transformation of the fluxon shape given by Eq.(18) (and probably by Eq.(14) as well ) takes 
place, this would result in a premature switching from the flux-flow branch and, possibly, in the existence of an extra 
metastable flux-flow branch in the IVC with the same limiting velocity, u = c\ , but with larger dissipation. 

The average DC voltage in J J I is V1/V01 = u/coi, where V01 = , and in J J 2 is zero. Therefore, we now plot 
current-velocity characteristics to represent IVC's. 

In Fig. 11, the single fluxon IVC's are shown for double stacks with equal damping coefficients, ai,2 = 0.05, and 
for Jcij J c \ = f i (solid diamonds); J c ij Jc\ — 0.5, (open circles); J c ij 1 J c \ — 2, (open squares). The rest of parameters 
are the same as in Fig. 1. In Fig. 11, symbols represent solutions obtained from ODE and subsequent solid lines 
show solution of the full PDE, Eq. (1), dashed gray line shows the IVC of an uncoupled single junction 1 and dotted 
line indicates the position of the lower Swihart velocity, c\ ~ 0.817coi. Insets in Fig. 11 show spatial distribution of 
sin(ipi) (solid lines), and sin(<fi2) (dashed lines), for maximum propagation velocities, for which ODE based numerical 
procedure converged: u max (J C 2/ J c i = 1) — 0.999ci, u max (J c2 / J c i = 0.5) ~ Q.97ei, u max (J c2 / J c i = 2) ~ 0.99ci. Solid 
lines in Fig. 11 show that PDE allow solutions propagating with even larger velocity, however those solutions are not 
of solitonic type, and will be discussed in the next section. From Fig. 11 it is seen that the IVC's for J c ij3c\ = 1 
and 0.5, exhibit velocity matching behavior at u — > c±. On the other hand, for J c ijJ c \ — 2, no velocity matching 
behavior is observed, and the velocity reaches c\ at a finite current. Such behavior is in agreement with the absence 
of contracted F\ component, as discussed above. From the insets in Fig. 11 it is seen, that the fluxon shape for the 
case of small damping does not differ much from the frictionless case, considered in the previous section. On the 
other hand, damping reduces the stability of the fluxon state at high propagation velocities, so that the maximum 
fluxon velocity for pure solitonic motion becomes less than c\. This is due to the finite bias current in the stack, 
which result in asymmetry of phase distribution in the stack. Such asymmetry is clearly seen from inset in Fig. 11 
for J C 2/Jci = 0.5. From Fig. 11 it is seen, that reduction of stability in the flux-flow state depends on parameters of 
the stack. E.g. for the case of double stack with nonidentical junctions, flux-flow state is more stable for fluxon in 
the weaker junction, J C 2/ Jci = 2, than for fluxon in the stronger junction, J c ij J c \ = 0.5. For the case J c il Jc\ = 0.5, 
switching of the second junction to the quasiparticle branch occurs first at I/I c i — 0.22. 

In real life of course junctions in the stack are not prenumerated and if junctions are not identical, the stable state 
will correspond to a fluxon placed in the weaker junction. Therefore whenever properties of junctions in the stack 
are considerably different, the stable dynamic state at u — c± would correspond to the existence of uncontracted F2 
component soliton in the weaker junction. Experimentally, steps at the velocity matching condition were observed 
for low-T c SJJ's jl8|,[l9]], on the other hand, for high-T c intrinsic SJJ's the flux-flow IVC's always have a finite slope 
ppf . We note, that for stacks with large number of junctions with thin electrodes, d <C A s , it is not necessary 
to have a variation in J C i to make the junctions different. In this case the middle junctions have a lower critical 
field, H c i 7 approximately half of that compared to the outmost junctions due to the fact that fluxon in the outmost 
junctions carries only half a flux quantum. In a double stack, considered here, a corresponding thing happens when 
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the junctions have different electrode thicknesses. The junction with thicker electrodes may have lower H c i even if J c 
of this junction is larger, see Ref. ||. In a sense, the criterion for transition from 'weak' to 'strong' junction is given 
by Eqs. (17,19). 

Although the stable state corresponds to the case when fluxon is placed in the weaker junction, the situation when 
the fluxon is placed in the stronger junction can also be achieved in experiment as shown in Ref. Jf9"| . Obviously such 
state can be achieved in annular SJJ's. In this case the fluxon can eventually be introduced in the stronger junction 
and if so it will stay there and can be accelerated by the bias current. 

C. Velocity above C\ : Cherenkov radiation. 

So far we have considered the case u < c\. For u > ci, coefficients before the second derivative of phase in Eq.(3) 
become negative. The equation can still be written in a sine-Gordon type form if we substitute tp(£) = ir + 0(f). 
Indeed, a 2ir soliton like solution for 0(f), moving with u > c%, does exist. However the energy of this state is 
decreasing with increasing velocity and therefore such state would correspond to unstable IVC branch with negative 
resistance. 

In Ref. H it was suggested that for u > c± the fluxon is a combination of a soliton with Josephson plasma waves. 
From Eq. (8) it is seen that Ai becomes imaginary for u > ci and the iq component transforms into a travelling 
Josephson plasma wave. The fluxon solution is then given by the Fi component soliton accompanied by Josephson 
plasma waves from the degenerate F± component. Recently, the existence of such type of solution was shown by 
numerical simulation in Ref. |2lf| and was interpreted as Cherenkov radiation in SJJ's, when fluxon velocity exceeds 
the phase velocity of electromagnetic waves. This solution is not of soliton type and can not be obtained from ODE. 
Therefore solution of full PDE, Eq.(l), is required. From Fig. 11 it is seen that for the case J C 2/J c i = 2, PDE allow 
solutions propagating with u > c\ . 

In Fig. 12 results of numerical simulations of PDE, Eq.(l), for the case J^ij J c \ = 2, u > c% are shown. Parameters 
of SJJ's are the same as in Fig. 11. Insets show spatial distributions of sin(ipi) (solid lines) and sin((f2) (dashed 
lines) for a) u/c\ ~ 1.015, I/I c \ = 0.15 and b) u/c\ ~ 1.155, I/I c i = 0.5, respectively. Simulations were done for 
annular SJJ's with L = IOOAji. From our simulations we observe that fluxon shape changes gradually as u exceeds 
ci. Therefore, there are no peculiarities at u = c\ in the IVC, see solid line in Fig. 11 for J c %l Jc\ = 2- Indeed, from 
inset a) in Fig. 12 it is seen that for u slightly above ci, fluxon shape in the left half-space is similar to that at u < ci, 
see the bottom inset in Fig. 11. However, small oscillations appear behind the fluxon (fluxon is propagating from right 
to left). As the velocity increases, both amplitude and wavelength of the oscillations increase, as illustrated in inset 
b) in Fig. 12. To clarify the physical origin of the oscillations, in Fig. 12 we have plotted the average wavelength of 
oscillations (circles) as a function of the absolute value of the Lorentz factor, y (u/ci) 2 — 1. The solid line in Fig. 12 
shows the absolute value, 27r |Ai|, given by Eq.(8) and describing small amplitude plasma waves from the degenerate 
F\ component Q| (the factor 2tt is due to different definition of the wavelength and the penetration depth). Excellent 
agreement between the wavelength of Cherenkov radiation and Josephson plasma wavelength from the degenerate F\ 
component is observed without any fitting, thus confirming the idea of Ref. j3[| that Cherenkov radiation is due to 
plasma wave generation from the degenerate Ft component. We would like to note, that | Ai | is not linear as a function 
of the Lorentz factor, although deviations from the linear dependence are small. For high propagation velocities, an 
uncertainty appears in determination of A. This is caused by the increase of the amplitude of oscillations, see inset b) 
in Fig. 12. Here oscillations are not exactly monochromatic but the wavelength slightly increase with the amplitude. 
The uncertainty in determination of A at high velocities is shown by error bars in Fig. 12. 

CONCLUSIONS 

In conclusion, we have shown that the shape of a single fluxon in double stacked Josephson junctions can be described 
by existence of two components given by Eqs. (6,7) and with characteristic lengths and velocities given by Eqs. (8,9). 
At velocities up to u ~ 0.98ci, the soliton shape is well described by the approximate double component solution, 
Eq.(10), for all studied junction parameters, as illustrated in Figs. 1-4, 7. In the very vicinity of the lower Swihart 
velocity the fluxon shape may undergo radical transformations. The final shape of the fluxon at u = ci strongly 
depends on parameters of the stack. From our numerical simulations we have found that the fluxon may remain 
double component, as shown in Fig.l, transform to a pure Fi component solution, see Fig. 3, or be a more complicated 
combination of F\ and Fi components, see Eqs. (14, 18) and Figs. 4,7,8. Those more complicated solutions do not, 
strictly speaking, represent the single fluxon state, but are combinations with fluxon-antifluxon pairs. However, even 
in this case these are always the components F\$ described by Eqs. (6,7) that constitute the solution. This implies 
that the components F\ i are real and may live their own life and appear in different combinations. Conditions for 
observing the single and the two-component solutions at u = c± were formulated and verified. We have shown that 
as the velocity approaches ci, the phase shift may become nonmonotonous. 

A prominent feature of a soliton moving at the velocity close to ci in SJJ's is the possible inversion of magnetic 
field £2(0), see Figs. 1,3, 4 c). Such behavior was predicted analytically in Ref. M. Here we confirmed the existence of 
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this phenomenon by numerical simulation. The inversion of magnetic held in junction 2 may lead to attractive fluxon 
interaction for fluxons in different junctions. Then the so-called 'in-phase' or 'bunched' state |1C(] with fluxons one on 
the top of the other in adjacent junctions may become favorable at high enough fluxon velocity. In experiment this 
would result in appearance of the extra flux-flow branch in the current-voltage characteristics with higher voltage, 
as shown by numerical simulations |l7| ] and observed in experiment on low-T c SJJ's |i~8| . In Rcf. JL4| it was shown 
that the bunched state can be stable at u > Ci, however, no mechanism for overcoming the mutual fluxon repulsion 
and transformation into the bunched state was suggested. The existence of the held inversion might be a criterion for 
the appearance of the bunched state in SJJ's. As we have shown, the sign inversion and a dip in -82(0) disappears 
when junction 2 becomes considerably stronger than junction 1 and transformation of the fluxon shape to a single F2 
component solution takes place, see Fig. 3 c). 

The shape of the flux-flow IVC's was analyzed for various parameters of SJJ's and it was shown that velocity 
matching behavior at u = Hi is observed when fluxon contains the contracted Fi component. Finally, Cherenkov 
radiation at u > c\ was shown to be due to generation of plasma waves from the degenerate F\ component in 
agreement with the prediction of Ref. Analytic expression for the wavelength of the Cherenkov radiation is 
derived. 

The work was supported by Swedish Superconductivity Consortium and in part by the Russian Foundation for 
Basic Research under Grant No. 96-02-19319. 
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FIG. 1. Profiles of a) phase differences v?i,2, b) the ratio sin(ipi) / sin(ip2) , and c) magnetic inductions Bi,2 of a single fluxon 
are shown for double SJJ's consisting of identical strongly coupled junctions and for different fluxon velocities, u/ci=0, 0.61, 
0.92, 0.98, 0.998, 0.9999 (from left to right curve). In Fig.l a) dotted lines show profiles obtained from the analytic double 
component solution Eq.(10). The rest of the curves represent results of numerical simulation. It is seen that the fluxon shape 
in this case is well described by Eq.(10) and consist of contracted and uncontracted components. The sign inversion of £2(0) 
at u ~ ci is clearly seen. 



FIG. 2. Fluxon shape at u = 0.98ci, for different critical currents, J C 2/ Jci~W', 2; 1; 0.5; 0.1, from left to right curve. 
The rest of the parameters of the stack and the way of presentation is the same as in Fig. la). It is seen that fluxon shape 
is in a good agreement with the analytical double component solution, Eq.(10), up to velocities very close to ci. A gradual 
transformation of the fluxon shape from the uncontracted F2 component solution to the contracted F\ component solution can 
be seen as J C 2/Jd decreases. 
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FIG. 3. Fluxon shape for the case when the fluxon is placed in the weaker junction, J C 2/Jd = 2, the rest of the parameters 
and the way of presentation are the same as in Fig.l. At velocities up to 0.98ci the shape of the fluxon is well described by the 
double component solution, Eq.(lO). At higher velocities transformation to the single F2 component solution, Eq. (6), with 
K2 = 2, takes place, as seen from Fig. 3 b). 

FIG. 4. Fluxon shape for the case when the fluxon is placed in the stronger junction, Jci/Jci = 0.5, the rest of the 
parameters and the way of presentation are the same as in Fig. 1. It is seen that at velocities up to 0.98ci the shape of the 
fluxon is well described by the double component solution, Eq.(10). At higher velocities the fluxon still has two contracted 
and uncontracted components Fi,2- The existence of the two fluxon components is clearly seen from Fig. 4 b). As in ci, 
sin(tpi) / sin(if2) — > K2 = 0.5 outside the center of the fluxon and sin(<pi) / sin(<p2) — ► —Ki = 1 in the center. However, 
transformation of the fluxon shape to that given by Eq.(14) takes place. 

FIG. 5. Phase distributions at u — 0.9999ci from Figs. 3a) and 4a). Solid and dashed curves represent ip\ and if 2, 
respectively, for the fluxon in the stronger junction J C 2j J c i = 0.5, while dashed-dotted and dotted curves represent ipi and <pz, 
respectively, for the fluxon in the weaker junction J^/Jd = 2. It is seen that <f2,i from Fig. 4a) merge with 931,2 from Fig. 3a) 
in the left half-space. Showing that the overall fluxon shape at u = ci for the case when the fluxon is in the stronger junction, 
Fig. 4, is given by Eq. (14). 

FIG. 6. Profiles of the fluxon moving with velocity very close to the lowest Swihart velocity, u = 0.9999ci, for different 
critical current densities J C 2/Jci = 10 — 0.1 increasing sequentially from the left to the right curve. The rest of the stack 
parameters and the way of presentation of data is the same as in Fig. 1. When the fluxon is placed in the weaker junction, 
the fluxon shape at the lowest Swihart velocity is described by the single F2 component. In Fig. 6 b) the dependence 
sin(ipi) I sin(ip2) — K2 = yy is clearly visible in the whole space region for J C 2/Jci > 1. When the fluxon is placed in the 
stronger junction, it has two components, Fi ; 2- From Fig. 6 b) it is seen that for the case J C 2/Jci < 1, the fluxon shape in the 
center is determined by the Fi component, sin(ipi) / sinfoz) = —m = 1, while outside the center the shape is given by the F2 
component, sin(tpi)/sin(ip2) = K2 = j^. 

FIG. 7. Fluxon shape for the case of nonidentical electrodes, A2 — 2.5Ai and J C 2/Jci = 0.5. The rest of the parameters 
and the way of presentation is the same as in Fig. 1. It is seen that at velocities up to 0.98ci the shape of the fluxon is 
well described by the double component solution, Eq.(10). From Fig. 7b) it is seen that outside the fluxon center the phase 
distribution is determined by the F2 component with K2 ~ 0.79 given by Eq.(13). However, at ti ~ 0.998ci a sudden switching 
to the other type of solutions occur. 

FIG. 8. Illustration of coexistence of three types of solution for the stack from Fig.7 for u ~ 0.998ci. a)fi,2 for the solution 
coming from low velocities, b) and c) show solutions given by Eq.(14)and Eq.(18), respectively, coming from higher velocities. 
Solid and dashed curves represent ip\ and (p2, respectively, obtained from numerical simulations, while dashed-dotted and 
dotted curves in Fig. 8 a) represent ip\ and <p2, respectively, given by the approximate double component solution, Eq.(10). 

FIG. 9. a) Soliton shape at u = 0.9999ci from Fig.7. Solid and dashed lines represent 931,2 obtained from numerical 
simulation, while dashed-dotted and dotted lines represent the analytic double component solution shifted by x — —xo and 
illustrating that the features at x = ±£0 correspond to the double component soliton in JJ 2. b) Inverted values of magnetic 
induction in the center of junction 1 versus the Lorentz factor. Symbols represent numerical simulations and solid line is an 
apparent linear fit. Clear Lorentz contraction of the central region at u = ci is seen, showing that it belongs to a pure Fi 
component. From Figs. 9 a,b) it is seen that the fluxon shape is given by Eq.(18). 

FIG. 10. Regions of existence of the two component (shaded area) and the single component solutions at u = c\ for ^ 2 ^ 2 = 1. 
Numbers show the number of components obtained numerically. Solid and dashed lines represent the conditions Eq.(l9) and 
Eq.(17), respectively. It is seen that the single component solution exists when both conditions are satisfied. Arrows indicate 
the cases considered in Figs. 1,3,4,7. 
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FIG. 11. Single fluxon IVC's are shown for double stack with equal damping, 0:1,2 = 0.05, and for J c2 / 'J c i —0.5, 1, and 2. 
The rest of parameters are the same as in Fig.l. Symbols and solid lines represent results obtained by solving ODE, Eq.(3) 
and PDE, Eq.(l), respectively. Dashed gray line shows the IVC of a single J J 1 and dotted line shows the position of the lower 
Swihart velocity, ci. Insets show spatial distribution of sin(ipi,2), for maximum propagation velocities, 

FIG. 12. The average wavelength of Cherenkov oscillations (circles) is shown as a function of the absolute value of the 
Lorentz factor. The solid line represents wavelength 27r|Ai|, for plasma waves from the degenerate Fi component, given by 
Eq.(8). Insets show spatial distributions of sin(ipi) (solid lines) and sin(<p2) (dashed lines) for a) w/ci ~ 1.015, I/I c i = 0.15 
and b) w/ci ~ 1.155, I/I c i = 0.5, respectively. Arrows indicate points in the main graph for which profiles were obtained. 
From Fig.12 it is seen that Cherenkov radiation at u > ci is caused by plasma waves from the degenerate F\ component. 
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